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O ' ABSTRACT 

\ Observations of redshifted 21cm radiation from neutral hydrogen (HI) at high redshifts is 

an important future probe of reionization. We consider the Multi-frequency Angular Power 
Spectrum (MAPS) to quantify the statistics of the HI signal as a joint function of the angular 
multipole / and frequency separation Av. The signal at two different frequencies is expected 
to decorrelate as Av is increased, and quantifying this is particularly important in deciding the 
frequency resolution for future HI observations. This is also expected to play a very crucial 
. role in extracting the signal from foregrounds as the signal is expected to decorrelate much 

' faster than the foregrounds (which are largely continuum sources) with increasing Av. 

\Q \ In this paper we develop formulae relating MAPS to different components of the three 

C^) ■ dimensional HI power spectrum taking into account HI peculiar velocities. We show that the 

flat-sky approximation provides a very good representation over the angular scales of interest, 
and a final expression which is very simple to calculate and interpret. We present results 
for z = 10 assuming a neutral hydrogen fraction of 0.6 considering two models for the HI 
distribution, namely, (i) DM: where the HI traces the dark matter and (ii) PR: where the effects 
of patchy reionization are incorporated through two parameters which are the bubble size and 
the clustering of the bubble centers relative to the dark matter (bias) respectively. We find that 
while the DM signal is largely featureless, the PR signal peaks at the angular scales of the 
individual bubbles where it is Poisson fluctuation dominated, and the signal is considerably 
enhanced for large bubble size. For most cases of interest at I ~ 100 the signal is uncorrected 
beyond Av ~ 1 MHz or even less, whereas this occurs around ~ 0.1 MHz at I ~ 10 3 . The 
53 ■ Av dependence also carries an imprint of the bubble size and the bias, and is expected to be 

an important probe of the reionization scenario. Finally we find that the I range 10 3 — 10 4 is 
optimum for separating out the cosmological HI signal from the foregrounds, while this will 
be extremely demanding at I < 100 where it is necessary to characterize the Av dependence 
of the foreground MAPS to an accuracy better than 1%. 

Key words: cosmology: theory, cosmology: diffuse radiation, cosmology: large-scale struc- 
ture of universe 
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1 INTRODUCTION 

One of the major challenges in modern cosmology is to understand 
the reionization history of the Universe. Numerous attempts have 
been made in this regard to constrain the evolution of neutral hy- 
drogen density (for recent reviews, see Barkana & Loeb 2001; Fan, 
Carilli & Keating 2006; Choudhury & Ferrara 2006a). Analyzes 
of quasar absorption spectra suggest that almost all of the hydro- 
gen has been ionized around redshift z < 6 (Becker et al. 2001; 
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Fan et al. 2002). A different set of constraints come from the Cos- 
mic Microwave Background (CMB) observations. The CMB pho- 
tons scatter with free electrons (produced because of reionization) 
which result in the suppression of the intrinsic anisotropies; at the 
same time a polarization signal is generated from this scattering. 
Measurements of the electron scattering optical depth from recent 
CMB observations like WMAP (Spergel et al. 2006; Page et al. 
2006) imply that reionization started before z ~ 10. It thus seems 
from these two results that reionization is an extended and complex 
process occurring over a redshift range 6-15 (Choudhury & Ferrara 
2006b; Alvarez et al. 2006). However, there exist limitations of us- 
ing these observations to study the details of reionization process 
which would be occurring at redshifts « 10. At present there is no 
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indication of quasars at such high redshifts; on the other hand, the 
electron scattering optical depth is sensitive only to the integral his- 
tory of reionization and it may not be useful to study the progress of 
reionization with redshift. In fact, it has been shown that the CMB 
polarization power spectrum is weakly dependent on the details of 
the reionization history (Kaplinghat et al. 2003; Hu & Holder 2003; 
Haiman & Holder 2003; Colombo 2004), though weak constraints 
could be obtained from upcoming experiments such as PLANCK 1 . 

Perhaps the most promising prospect of studying reionization 
at various stages of its progress is through future 21 cm observa- 
tions(for recent review see Furlanetto et al. (2006)). This basically 
involves observing the redshifted 21 cm line from the neutral hy- 
drogen (HI) at high redshifts. The advantage of these experiments 
lie in the fact that one can track neutral hydrogen fraction at any 
desired redshift by appropriately tuning the observation frequency. 
Tomography of the HI distribution using observations of the red- 
shifted 21 cm line is one of the most promising tool to study reion- 
ization. 

The possibility of observing 21 cm emission from the cos- 
mological structure formation was first recognized by Sunyaev and 
Zeldovich (1972) and later studied by Hogan & Rees (1979), Scott 
& Rees (1990) and Madau, Meiksin & Rees ((1997)) considering 
both emission and absorption against the CMB. More recently, the 
effect of heating of the HI gas and its reionization on 21 cm sig- 
nal has drawn great deal of attention and has been studied in detail 
(Gnedin & Ostriker 1997; Shaver et al. 1999; Tozzi et al. 2000; 
Iliev et al. 2002, 2003; Ciardi & Madau 2003; Furlanetto, Sokasian 
& Hernquist 2004; Miralda-Escude 2003; Chen & Miralda-Escude 
2004; Cooray & Furlanetto 2005; Cooray 2005; Mcquinn et al. 
2005; Sethi 2005; Salvaterra et al. 2005; Carilli 2006). 

It is currently perceived that a statistical analysis of the fluctu- 
ations in the redshifted 21 cm signal which is present as a minute 
component of the background radiation in all low frequency ra- 
dio observations holds the greatest potential. This approach has 
been considered in the context of lower redshifts (Bharadwaj, Nath 
& Sethi 2001; Bharadwaj & Sethi 2001; Bharadwaj & Pandey 
2003; Bharadwaj & Srikant 2004). A similar formalism can also 
be applied at high redshifts to probe reionization and also the pre- 
reionization era (Zaldarriaga, Furlanetto & Hernquist 2004; Furlan- 
etto, Zaldarriaga & Hernquist 2004b; Morales & Hewitt 2004; 
Bharadwaj & Ali 2004, 2005; Bharadwaj & Pandey 2005; Ali, 
Bharadwaj and Pandey 2005, 2006; Loeb & Zaldarriaga 2004; He 
et al. 2004). One of the main aims of epoch of reionization 21 cm 
studies is to probe the size, spatial distribution and evolution of 
the ionized regions which has been addressed by several groups 
(Furlanetto, McQuinn & Hernquist 2006; Furlanetto, Zaldarriaga 
& Hernquist 2004a; Wyithe & Loeb 2004; Mellema et al. 2006; 
Iliev et al. 2005). 

On the observational front, several initiatives are currently un- 
derway. The Giant Metre- Wave Radio Telescope (GMRT 2 ; Swarup 
et al. 1991) is already functioning at several bands in frequency 
range 150-1420 MHz and can potentially detect the 21 cm sig- 
nal at high redshifts. In addition, upcoming low-frequency exper- 
iments such as LOw Frequency ARray (LOFAR 3 ), Square Kilo- 
meter Array (SKA 4 ), PrimevAl Structure Telescope (PAST 5 ) have 
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now raised the possibility to detect 21 cm signal from very high 
redshift and thus motivated a detailed study of the expected 21 cm 
background from neutral hydrogen at high redshifts. 

Although the redshifted 21 cm line can provide enormous 
amounts of information, its detection is going to be a huge chal- 
lenge. The signal is expected to be highly contaminated by fore- 
ground radio emission. Potential sources for these foregrounds in- 
clude synchrotron and free-free emission from or Galaxy and exter- 
nal galaxies, low-frequency radio point sources and free emission 
from electrons in the intergalactic medium (Shaver et al. 1999; Di- 
Matteo. et al. 2002, 2004; Oh 1999; Cooray & Furlanetto 2004). 
Since these foregrounds are much stronger than the 21 cm signal 
arising from HI, it has been suggested that it may be impossible to 
study reionization via such observations (Oh & Mack 2003). 

However, there have been various proposals for tackling 
such foregrounds, the most promising being the usage of multi- 
frequency observations. It has been proposed that multi-frequency 
analysis of the radio signal can be useful in separating out the fore- 
ground (see e.g. Shaver et al. 1999; DiMatteo. et al. 2002; Gnedin 
& Shaver 2004; DiMatteo. et al. 2004; Santos, Cooray & Knox 
2005). In fact this has also been noted by most of the authors who 
have studied the statistical signal from redshifted HI and have been 
referred to earlier. The cross correlation between the HI signal at 
different frequencies is expected to decay rapidly as the frequency 
separation is increased (Bharadwaj & Ali 2005; Santos, Cooray & 
Knox 2005) while the foregrounds are expected to have a contin- 
uum spectrum and hence should be correlated across frequencies. 
This property of different foregrounds contaminants can in princi- 
ple be used to remove them from expected 21 signal. 

In this paper, we develop the formalism to calculate the multi- 
frequency angular power spectrum (hereafter MAPS) which can be 
used to Analise the 21 cm signal from HI both in emission and 
absorption against the CMB. We restrict our attention to HI emis- 
sion which is the situation of interest for the epoch of reionization. 
In our formalism, we consider the effect of redshift-space distor- 
tions which has been ignored in many of earlier works. As noted 
by Bharadwaj & Ali (2004), this is an important effect and can en- 
hance the mean signal by 50% or more and the effect is expected 
to be most pronounced in the multi-frequency analysis. We next 
use the flat sky approximation to develop a much simpler expres- 
sion of MAPS which is much easier to calculate and interpret than 
the angular power spectrum written in terms of the spherical Bessel 
functions. We adopt a simple model for the HI distribution (Bharad- 
waj & Ali, 2005) which incorporates patchy reionization and use it 
to predict the expected signal and study its multi-frequency prop- 
erties. The model allows us to vary properties like the size of the 
ionized regions and their bias relative to the dark matter. We use 
MAPS to analyze the imprint of these features on the HI signal and 
discuss their implication for future HI observations. 

As noted earlier, the HI signal at two different frequencies 
separated by Av is expected to become uncorrelated as Av is in- 
creases. As noted in Bharadwaj & Ali (2005), the value of Au 
beyond which the signal ceases to be correlated depends on the 
angular scales being observed and it is < 1 MHz in most situa- 
tions of interest. A prior estimate of the multi-frequency behavior 
is extremely important when planning HI observations. The width 
of the individual frequency channels sets the frequency resolution 
over which the signal is averaged. This should be chosen suffi- 
ciently small so that the signal remains correlated over the channel 
width. Choosing a frequency channel which is too wide would end 
up averaging uncorrelated HI signal which would wash out vari- 
ous important features in the signal, and also lead to a degradation 
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in the signal to noise ratio. In this context we also note that an 
earlier work (Santos, Cooray & Knox, 2005) assumed individual 
frequency channels 1 MHz wide and smoothed the signal with this 
before performing the multi-frequency analysis. This, as we have 
already noted and shall study in detail in this paper, is considerably 
larger than the Av where the signal is uncorrected and hence is 
not the optimal strategy for the analysis. We avoid such a pitfall by 
not incorporating the finite frequency resolution of any realistic HI 
observations. It is assumed that the analysis of this paper be used 
to determine the optimal frequency channel width for future HI ob- 
servations. Further, it is quite straightforward to introduce a finite 
frequency window into our result through a convolution. 

The outline of this paper is as follows. In Section 2 we present 
the theoretical formalism for calculating MAPS of the expected 21 
cm signal considering the effect of HI peculiar velocity. The calcu- 
lation in the full-sky and the fiat-sky approximation are both pre- 
sented with the details being given in separate Appendices. Section 
2.1 defines various components of the HI power spectrum and Sec- 
tion 2.4 presents to models for the HI distribution. We use these 
models when making predictions for the expected HI signal. We 
present our results in Section 3 and also summarize our findings. 
In Section 4. we discuss the implications for extracting the signal 
from the foregrounds. 



2 THEORETICAL FORMALISM 
2.1 The HI power spectrum 

The aim of this section is to set up the notation and calculate the 
angular correlation Ci for the 21cm brightness temperature fluctu- 
ations. It is now well known (e.g.. Bharadwaj & Ali 2005) that the 
excess brightness temperature observed at a frequency v along a 
direction n is given by 



T(u,n) =T(z) m i (z,nrv) 



(1) 



where the frequency of observation is related to the redshift by v = 
1420/(1 + z) MHz. The quantity 
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is the comoving distance and 
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(3) 



where Y w 0.24 is the helium mass fraction and all other sym- 
bols have usual meaning. In the above relation, it has been assumed 
that the Hubble parameter H(z) « HoQU. 2 (l + z) 3 ^ 2 , which is a 
good approximation for most cosmological models at z > 3. The 
quantity r/m is known as the "21 cm radiation efficiency in redshift 
space" (Bharadwaj & Ali 2005) and can be written in terms of the 
mean neutral hydrogen fraction xm and the fluctuation in neutral 
hydrogen density field 8m as 



rjm(z, nr„) 



xm(z)[l + Sm(z,nr v )] ( 1 
(1 + z) dv(z, nr„) 



H(z) dr v 



(4) 



where T 7 and T s are the temperature of the CMB and the spin 
temperature of the gas respectively. The term in the square bracket 
arises from the coherent components of the HI peculiar velocities. 
At this stage, it is useful to make a set of assumptions which 



will simplify our analysis: (i) We assume that T s 3> T 7 , which 
corresponds to the scenario where the intergalactic gas is heated 
well above the CMB temperature during the reionization process. 
This assumption is expected to be valid throughout the IGM soon 
after the formation of first sources of radiation as a substantial back- 
ground of X-rays from supernovae, star-formation, etc. can heat the 
IGM quickly, (ii) We assume that the HI peculiar velocity field is 
determined by the dark matter fluctuations, which is reasonable as 
the peculiar velocities mostly trace the dark matter potential wells. 
We then have 



r)m(z, nrv) = 



/ 



d 3 k 
(2^)3 



ifcrv (k-n) 



Vm (z,k) 



where 



fjm (z, k) = xm(z) [a H i(z, k) + (k ■ n) 2 A(z, k) 



(5) 



(6) 



and Ahi(z, k) and A(z, k) are the Fourier transform of the fluc- 
tuations in the HI and the dark matter densities respectively. Note 
that /(fl m ), which relates peculiar velocities to the dark matter, 
has been assumed to have a value f(Sl m ) = 1 which is reasonable 
at the high z of interest here. 

For future use, we define the relevant three dimensional (3D) 
power spectra 

<A(z,k)A*(z,k')> = (2nfS D (k-k')P(z,k) 
<A H i(z,k)AHi(z,k')> = (27r) 3 fo(k-k')P A 2 i ( Z ,fc) 

<A( z ,k)AHi(z,k')} = (27r) 3 fo(k-k')PA HI (z,fc) (7) 

where P(z, k) and P A 2 (z, k) are the power spectra of the fluc- 
tuations in the dark matter and the HI densities respectively, while 
•Pahi (z, k) is the cross-correlation between the two. 



2.2 The multi-frequency angular power spectrum (MAPS) 

The multi-frequency angular power spectrum of 21 cm brightness 
temperature fluctuations at two different frequencies v\ and h>2 is 
defined as 



Ci(vx,v<i) EE (ai m {vi) a* m (z/ 2 )) ■ 



(8) 



In our entire analysis v\ and vi are assumed to differ by only a 
small amount Av <C v\, and it is convenient to introduce the nota- 
tion 



d(Ais) = C,(i/,i/ + Ai/) 



(9) 



where we do not explicitly show the frequency v whose value will 
be clear from the context. Further, wherever possible, we shall not 
explicitly show the z dependence of various quantities like T, xm, 
P(k) etc., and it is to be understood that these are to be evaluated 
at the appropriate redshift determined by v. 

The spherical harmonic moment of T(y, n) are defined as 



ai m {v) 



J dfil7 m (n)2>,n) 



kr u (k-n) 



(10) 



Putting the expression (6) for j^hi (k) in the above equation, one 
can explicitly calculate the MAPS in terms of the three dimensional 
power spectra defined earlier. We give the details of the calculation 
in Appendix A and present only the final expression for the angular 
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power spectrum at a frequency v 

Ci(Av) = HI / fc 2 dfc \ji(kr„)ji(k 



»V 2 )P A 2 [i (fe) 



- {ji{kr v )j'i \kr V2 ) + ji(kr V2 )j" (kr v )} P Alil (k) 
+ j'/ikr^j'/ikr^Pik)] 



(11) 



Here = -^riji{x) and we have used the notation r„ 2 

r v + r£, A;/ with 



dr v 



c (1 + z) 2 
~vo H{z) 



(12) 



Note that equation (11) predicts C; (Av) from the cosmological 21 
cm HI signal to be real. 

With increasing Av, we expect the two spherical Bessel func- 
tions ji (kr vi ) and ji (kr V2 ) to oscillate out of phase. As a con- 
sequence the value of Ci(Av) is expected to fall increasing Av. 
We quantify this through a dimensionless frequency decorrelation 
function defined as the ratio 



m(Av) 



Ci(Q) 



(13) 



For a fixed multipole /, this fall in this function with increasing 
Av essentially measures how quickly features at the angular scale 
9 ~ it /I in the 21 cm HI maps at two different frequencies become 
uncorrelated. Note that < |k;(A^)| < 1. 



2.3 Flat-sky approximation 

Radio interferometers have a finite field of view which is deter- 
mined by the parameters of the individual elements in the array. 
For example, at 150 MHz this is around 3° for the GMRT. In most 
cases of interest it suffices to consider only small angular scales 
which correspond to I >• 1. For the currently favored set of flat 
ACDM models, a comoving length scale R at redshift z > 5 would 
roughly correspond to a multipole 



I « 3 x 10 4 1 - 



1.1 



R 



(14) 



VI + zj \h- 1 Mpc / 

Thus, for length scales of R < lOO/t -1 Mpc at z ^ 10, one would 
be interested in multipoles I > 200. For such high values of I one 
can work in the flat-sky approximation. 

A small portion of the sky can be well approximated by a 
plane. The unit vector n towards the direction of observation can 
be decomposed as 



n = m + 0;m-0 = O; I 9 |« 1 



(15) 



where m is a vector towards the center of the field of view and 8 is 
a two-dimensional vector in the plane of the sky. It is then natural 
to define the two-dimensional Fourier transform of T(v, n) in the 
flat-sky as 



f(v,XJ) 



(16) 



where U, which corresponds to an inverse angular scale, is the 
Fourier space counterpart of 6. Using equations (1) and (5), and the 
fact that for the flat-sky we can approximate k • n « k • m = k\\ 
we have 



7>,U) = 



T 



2tt r 



if 



dfc„ e ifc " r - 7?Hi(fc||m + 27rU/rv) . (17) 



It is useful to introduce the rjm power spectrum Phi defined as 

(j?Hi(k)j7Hi(k')) = (27r) 3 4(k - k')Pm(k) . (18) 

This is related to the other three power spectra introduced earlier 
through 

Pm(k) = xl^lP^ik) +2 M 2 P Am (fc) +fJ. 4 P(k)] (19) 

where /x = m • k = fcy/fc. Note that the anisotropy of Pni(k) i.e., 
its /u-dependence arises from the peculiar velocities. 

The quantities calculated in the flat-sky approximation can be 
expressed in terms of their all-sky counterparts. The correspon- 
dence between the all-sky angular power spectra and its flat-sky 
approximation is given by 

(TVi,U)T> 2 ,U')> =CW iyi,V2) 4 2) (U-U') (20) 

where <5^(U — U') is the two-dimensional Dirac-delta function. 
The details of the above calculation are presented in Appendix B. 
Thus allows us to estimate the angular power spectrum Ci under 
the flat-sky approximation which has a much simpler expression 
(Bharadwaj & Ali 2005) 



C*? at (A 



rpl f- OC 

v) = —\ dfc|, cos(fc||7vAz/)PHi(k) (21) 



where the vector k has magnitude k = ^Jk'j + I 2 /r% ie. k has 

components fcii and l/r v along the line of sight and in the plane 
of the sky respectively. It is clear that the angular power spectrum 
Ci(Av) is calculated by summing over all Fourier modes k whose 
projection in the plane of the sky has a magnitude l/r v . We also 
see that C'i is determined by the power spectra only for modes k > 
J/r„. 

The flat-sky angular power spectrum Cf at (0) is essentially 
the 2D power spectrum of the HI distribution on a plane at the dis- 
tance r v from the observer, and for Av = equation (21) is just 
the relation between the 2D power spectrum and its 3D counterpart 
(Peacock 1999). For Av 7^ it is the cross-correlation of the 2D 
Fourier components of the HI distribution on two different planes, 
one at r v and another at r v +Av Any 2D Fourier mode is calculated 
from its full 3D counterparts by projecting the 3D modes onto the 
plane where the 2D Fourier mode is being evaluated. The same set 
of 3D modes contribute with different phases when they are pro- 
jected onto two different planes. This gives rise to cos(fc||r^A^) 
in equation (21) when the same 2D mode on two different planes 
are cross-correlated and this in turn causes the decorrelation of 
Cf at (A*/) with increasing Av. 

Testing the range of I over which the flat-sky approximation is 
valid, we find that for the typical HI power spectra Cf at ( Av) is in 
agreement with the full-sky Ci (Av) calculated using equation (11) 
at a level better than 1 per cent for angular modes / > 10. Since 
the integral in equation (21) is much simpler to compute, and more 
straightforward to interpret, we use the flat-sky approximation of 
Ci for our calculations in the rest of this paper. 

Note that equation (21) is very similar to the expression for 
the visibility correlations [equation (16) of Bharadwaj & Ali 2005] 
expected in radio interferometric observations of redshifted 21 cm 
HI emission. The two relations differ only in a proportionality fac- 
tor which incorporates the parameters of the telescope being used 
for the observation. This reflects the close relation between the vis- 
ibility correlations, which are the directly measurable quantities in 
radio interferometry, and the Cis considered here. 
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2.4 Modeling the HI distribution 

The crucial quantities in calculating the angular correlation func- 
tion are the three dimensional power spectra P(k), PA UI (k) and 
P A 2 (k). The form of the dark matter power spectrum P(k) is rel- 
atively well-established, particularly within the linear theory. We 
shall be using the standard expression given by Bunn & White 
(1996). 

The power spectrum of HI density fluctuations P A 2 (k) and 
its cross-correlation with the dark matter fluctuations Pa hi (k) are 
both largely unknown, and determining these is one of the most im- 
portant aims of the future redshifted 21 cm observations. We adopt 
two simple models with a few parameters which capture the salient 
features of the HI distribution. 

The first model, which we shall denote as DM, assumes ho- 
mogeneous reionization where the HI traces the dark matter , i.e., 
Ahi = A. In this model we have 

Pt&W = P± m (k) = P(k) (22) 

which we use in equation (1 1) to calculate Ci (Ai/). Alternately, we 
have 

Pm(fc) = x 2 HI (l + fJ 2 ) 2 P(k) (23) 

which we can use in equation (21) to calculate Ci(Av) in the flat- 
sky approximation. This model has only one free parameter namely 
the mean neutral fraction xm- 

The second model, denoted as PR, incorporates patchy reion- 
ization. It is assumed that reionization occurs through the growth 
of completely ionized regions (bubbles) in the hydrogen distribu- 
tion. The bubbles are assumed to be spheres, all with the same co- 
moving radius R, their centers tracing the dark matter distribution 
with a possible bias b c . While in reality there will be a spread in 
the shapes and sizes of the ionized patches, we can consider 7? as 
being the characteristic size at any particular epoch. The distribu- 
tion of the centers of the ionized regions basically incorporates the 
fact that the ionizing sources are expected to reside at the peaks 
of the dark matter density distribution and these are expected to 
be strongly clustered. For non-overlapping spheres the fraction of 
ionized volume is given by 

-, - 4ttP 3 _ 
xmi = l-xm = — - — ram (24) 
o 

where nm is the mean comoving number density of ionized 
spheres and we use TZ to denote the ratio TZ = xmi/xm- This 
model has been discussed in detail in Bharadwaj and Ali (2005), 
and we have 

Am(k) = [1 - b c TZW(kR)} A(k) - TZW(kR)A P (k) . (25) 

The HI fluctuation is a sum of two parts, one which is corre- 
lated with the dark matter distribution and an uncorrelated Pois- 
son fluctuation Ap. The latter arises from the discrete nature of 
the HII regions and has a power spectrum ra^j 1 . Also, W(y) = 
(3/j/ 3 )[siny — j/cosy] is the spherical top hat window function 
arising from the Fourier transform of the spherical bubbles. This 
gives 

P A 2 (k) = [1 - TZb c W{kR)] 2 P(k) + ^ nW ^ kR )f (26 ) 
HI ram 

and 

Pahi (k) = [1 - TZb c W{kR)] P{k) (27) 
which we use in equation (1 1) to calculate Ci {Av). Alternately, we 



have the HI power spectrum (Bharadwaj and Ali, 2005) 

Pm(k) = x m f [1 - TZb c W(kR) + ^ 2 ] 2 P(k) + } 

(28) 

which we can use in equation (21) to calculate Ci(Av) in the flat- 
sky approximation. 

Our analysis assumes non-overlapping spheres and hence it is 
valid only when a small fraction of the HI is ionized and the bias 
is not very large. As a consequence we restrict these parameters to 
the range xm > 0.5 and b c < 1.5. We note that in the early stages 
of reionization (ie. ami <C 1) equation (28) matches the HI power 
spectrum calculated by Wang & Hu (2005), though their method of 
arriving at the final result is somewhat different and is quite a bit 
more involved. 

For large values of the bubble size (R > 8 Mpc) the power 
spectrum is essentially determined by the Poisson fluctuation term. 
At length-scales larger than the bubble size (k < n/R) we have 
W 2 (kR) « 1, and hence the power spectrum is practically con- 
stant. Around scales corresponding to the characteristic bubble size 
k w n/R, the window function W(kR) starts decreasing which 
introduces a prominent drop in Pm(k). For smaller length-scales 
(k > n/R), the power spectrum shows oscillations arising from the 
window function W(kR). At sufficiently large k the power spec- 
trum Pm(k) is dominated by the dark matter fluctuations and it 
approaches the DM model. We should mention here that the oscil- 
lations in Pm(k) are a consequence of the ionized bubbles being 
spheres, all of the same size. In reality there will be a spread in the 
bubble shapes and sizes, and it is quite likely that such oscillation 
will be washed out (Wang & Hu 2005) but we expect the other fea- 
tures of the PR model discussed above to hold if the characteristic 
bubble size is large (R > 8 Mpc). 

For smaller values of R, the power spectrum Pm(k) could 
be dominated either by the term containing the dark matter power 
spectrum P(k) or by the Poisson fluctuation term, depending on 
the value of lZb c — 2b c /3. 

In addition to the effects considered above, the random mo- 
tions within clusters could significantly modify the signal by elon- 
gating the HI clustering pattern along the line of sight [the Finger 
of God (FoG) effect]. We have incorporated this effect by multi- 
plying the power spectrum Phi (k) with an extra Lorentzian term 
(l + fcjjo-fVa 2 ^)- 1 (Sheth 1996; Bollinger et al. 1996) where cr P 
is the one dimensional pair-velocity dispersion in relative galaxy 
velocities. 

The cosmological parameters used throughout this paper are 
those determined as the best-fit values by WMAP 3-year data 
release, i.e., Q m = 0.23, n b h 2 = 0.022, n s = 0.96, ft = 
0.74, a s = 0.76 (Spergel et al. 2006). Further, without any loss 
of generality, we have restricted our analysis to a single redshift 
z = 10 which corresponds to a frequency v — 129 MHz, and have 
assumed Shi = 0.6 (implying TZ = 2/3) which is consistent with 
currently favored reionization models. 

The effects of patchy reionization are expected to be most 
prominent at around z w 10 in currently favored reionization mod- 
els. At higher redshifts, the reionization is in its preliminary stages 
(xmi <S 1) and the characteristic bubble size R is quite small. We 
expect the DM model to hold at this epoch. At lower redshifts the 
HI signal is expected to be drastically diminished because most 
of the hydrogen would be ionized. Given this, it is optimum to 
study the HI signal properties at some intermediate redshift where 
xhii ~ 0.5 and R is reasonably large. For the currently favored 
reionization scenarios, it seems that these properties are satisfied at 
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Figure 1. The angular power spectrum of HI brightness temperature fluc- 
tuations for different models of the HI distribution assuming b c = 1 and 
CTp = 0. 

z ~ 10 (Choudhury & Ferrara 2006b), which we shall be studying 
in the rest of this paper. 



3 RESULTS 

We first consider the angular power spectrum Ci(Av) at Av = 
for which the results are shown in Figure 1. As discussed earlier, 
Ci (0) is essentially the 2D power spectrum of HI fluctuations eval- 
uated at the 2D Fourier mode l/r v ~ I x 10 - Mpc -1 . The results 
for the DM model serve as the fiducial case against which we com- 
pare different possibilities for patchy reionization. 

For large bubble size (R > 8 Mpc) the HI signal is dominated 
by Poisson fluctuations and it is well described through 

+ l)Ci oc y/^Tif — I (29) 
r v 

on scales larger than the bubble. At these angular scales the HI 
signal is substantially enhanced compared to the DM model. For 
smaller bubble size, the large angle signal is sensitive to the bias 
b c . The signal is very similar to the DM model for b c — and it is 
suppressed for higher bias. In all cases (large or small bubble size), 
the signal is Poisson fluctuation dominated on scales comparable to 
the bubble size and it peaks at I « nr^/R, with no dependence on 
b c . The HI signal traces the dark matter on scales which are much 
smaller than the bubble size. 

We next consider the behavior of Ki(Av), the frequency 
decorrelation function shown in Figure 2. For the DM model (upper 
left panel) where the HI fluctuations trace the dark matter we find 
that the frequency difference Av over which the HI signal remains 
correlated reduces monotonically with increasing I. For example, 
while for I = 100 K t (Av) falls to ~ 0.5 at Av ~ 500 KHz, it 
occurs much faster {Av ~ 10 KHz) for I — 10 . Beyond the first 
zero crossing ki (Av) becomes negative (anti-correlation) and ex- 
hibits a few highly damped oscillations very close to zero. These 
oscillations arise from the cos term in equation (21). The change 
in the behavior of ki (Av) for the DM model arising from the FoG 
effect is also shown in the same panel. Wang & Hu (2005) have 
proposed that a v is expected to have a value ~ 30 km/s at z ~ 8; 
in view of this, we show results for a p = 20 and 40 km/s. We 
find that there is a discernible change at I > 10 4 , and the FoG ef- 
fect causes the signal to remain correlated for a larger value of Av. 
For rj p = 20 km/s, the change is at most 15% for I = 10 4 and 
around 100% at I = 10 s . Though we have not shown it explicitly, 
we expect similar changes due to FoG effect in the PR model also. 



The patchy reionization model shows distinct departures from 
the DM model in the behavior of ki(Av). This reflects the im- 
print of the bubble size and the bias on the Av dependence. For 
R = 2 Mpc and b c = 1 (upper right panel) the large I (I > 1000, 
comparable to bubble size) behavior is dominated by the Poisson 
fluctuation of the individual bubbles which makes ni(Av) quite 
distinct from the DM model. Notice that for I = 10 3 , ki (Av) falls 
faster than the DM model whereas for I — 10 4 it falls slower than 
the DM model causing the / = 10 3 and 10 4 curves to nearly over- 
lap. The oscillations seen in Ci as a function of I in Figure 1 are also 
seen in the Av dependence of Ki (Av) at large I (10 5 ). The behav- 
ior at I = 10 2 is a combination of the dark matter and the ionized 
bubbles, and is sensitive to & c . For b c = 1, the initial decrease in 
ki (Av) is much steeper than the DM model with a sudden break 
after which the curve flattens. Figure 3 shows the b c dependence for 
/ = 10 2 and 10 3 . The bias dependence is weak for / = 10 3 where 
the Poisson fluctuations begin to dominate. For I — 10 2 , changing 
b c has a significant affect only near the break in ki(Av) leaving 
much of the curve unaffected. For a smaller bubble size we ex- 
pect a behavior similar to R = 2Mpc, with the b c dependence be- 
ing somewhat more pronounced and the Poisson dominated regime 
starting from a larger value of I. 

For large bubble size (R > 8Mpc) the large angle HI signal 
(I < ivr^/R) is entirely determined by the Poisson fluctuations 
where the signal is independent of I. This is most clearly seen for 
R = 8Mpc where the ki (Av) curves for I — 10 2 and / = 10 3 are 
identical. For both R — 8Mpc and 32 Mpc the large I behavior of 
ki(Av) approaches that of the DM model. 

In the final part we quantify the frequency difference Av 
across which the HI signal at two different frequencies remain cor- 
related. To be more precise, we study the behavior of Av^m which 
is defined such that ki(Av 1 / 2 ) ~ 1/2 ie. the correlation falls to 
50% of its peak value at Av = 0. We study this for different an- 
gular scales (different /) for the various models of HI distribution 
considered here. The main aim of this exercise is to determine the 
frequency resolution that would be required to study the HI fluc- 
tuations on a given angular scale /. Optimally one would like to 
use a frequency resolution smaller than Avx/2- A wider frequency 
channel would combine different uncorrelated signals whereby the 
signal would cancel out. Further, combining such signals would not 
lead to an improvement in the signal to noise ratio. Thus it would 
be fruitful to combine the signal at two different frequencies only 
as long as they are correlated and not beyond, and we use Av\ii to 
estimate this. The plot of Av x / 2 vs I for the different HI models is 
shown in Figure 4. 

We find that for the DM model Avxn falls monotonically 
with I and the relation is well approximated by a power law 

Av 1/2 = 0.2 MHz x f— J (30) 

which essentially says that A^ 1( / 2 ~ 0.66 MHz on 1° angular 
scales, Avi/2 ~ 0.04 MHz on V angular scales and Avx/ 2 ~ 
2 KHz on 1" angular scales. 

For the PR model, as discussed earlier, K; (Av) is I indepen- 
dent on angular scales larger than the bubble size I < nr v /R. As 
a consequence Av x /2 also is independent of I and it depends only 
on the bubble size R. This can be well approximated by 

Av 1/2 w 0.04 MHz fr-P-] (31) 
\Mpc / 

which given a large value 1.3 MHz for R — 32 Mpc while it falls 
below the DM model for 0.08 MHz for R = 2 Mpc. The large 
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Figure 2. The frequency decorrelation function /t;(Ai/) [defined in equation (13)] at £ = 10 2 , 10 3 , 10 4 , 10 5 . Results are shown for the DM model and the 
PR model with b c = 1 and the R values shown in the figure. For the DM model, we show results incorporating the FoG effect using ap = 20 and 40 km/a. 
For each I value re; ( Au) decreases faster for ap = and slowest for ap = 40 km / s. There is a significant change due to the FoG effect only at / > 10 4 . 
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Figure 4. This shows Av-^ji vs. / for the DM model and the PR model with b c = 1 for the R values shown in the figure. 



I behavior of Az/1/2 approaches the DM model though there are 
oscillations which persist even at large I. 

We note that our findings are consistent with the earlier find- 
ings of Bharadwaj & Ali (2005) whereas they significantly differ- 
ent from the results of Santos, Cooray & Knox (2005) who assume 
frequency channels of 1 MHz which is too large. 



4 IMPLICATIONS FOR SEPARATING SIGNAL FROM 
FOREGROUNDS 

Astrophysical foregrounds are expected to be several order of mag- 
nitude stronger than the 21 cm signal. The MAPS foreground con- 
tribution at a frequency v can be parametrized as (Santos, Cooray 
& Knox, 2005) 



Ci(Av) = A 



(?)" 



v + Ai 



Ii(Av) 



(32) 



where Vf — 130 MHz and a is the mean spectral index. The 
actual spectral index varies with line of sight across the sky and 
this causes the foreground contribution to decorrelate with in- 
creasing frequency separation Av which is quantified through the 
foreground frequency decorrelation function Ii(Av) (Zaldarriaga, 
Furlanetto & Hernquist, 2004) which has been modeled as 



Ii(Av) = exp 



logi 



v ' 



(33) 



We consider the two most dominant foreground components 
namely extragalactic point sources and the diffuse synchrotron ra- 
diation from our own galaxy. Point sources above a flux level S cu t 
can be identified in high-resolution images and removed. We as- 
sume Scut = O.lmJy and adopt the parameter values from Ta- 
ble 1 of Santos, Cooray & Knox (2005) for A, a, (3 and £. Figure 
5 shows the expected C*;(0) for the signal and foregrounds. The 
galactic synchrotron radiation dominates at large angular scales 
I < 10, 000 while the extragalactic point sources dominate at small 
angular scales. For all values of I, the foregrounds are at least two 
orders of magnitude larger than the signal. 

The foregrounds have a continuum spectra, and the contribu- 
tions at a frequency separation Av are expected to be highly cor- 
related. For Av = 1 MHz, the foreground decorrelation function 
7; ( Av) falls by only 2 x 10~ 6 for the galactic synchrotron radiation 
and by 3 x 10~ 5 for the point sources. In contrast, the HI decorrela- 
tion function ki(Av) is nearly constant at very small Av and then 
has a sharp drop well within 1 MHz, and is largely uncorrelated be- 
yond. This holds the promise of allowing the signal to be separated 



from the foregrounds. A possible strategy is to cross-correlate dif- 
ferent frequency channels of the full data which has both signal and 
foregrounds, and to use the distinctly different Av dependence to 
separate the signal from the foregrounds (Zaldarriaga, Furlanetto & 
Hernquist, 2004). An alternate approach is to subtract a best fit con- 
tinuum spectra along each line of sight (Wang et al., 2006) and then 
determine the power spectrum. This is expected to be an effective 
foreground subtraction method in data with very low noise levels. 
We consider the former approach here, and discuss the implications 
of our results. 

MAPS characterizes the joint / and Av dependence which 
is expected to be different for the signal and the foregrounds. 
For a fixed I, it will be possible to separate the two with rela- 
tive ease at a frequency separation Av if the decrement in the 
signal C'i (0)[1 — ni(Av)] is more than that of the foregrounds 
Ci (0)[1 — Ii(Av)]. Note that because the foregrounds are much 
stronger than the HI signal, a very small decorrelation of the fore- 
ground contribution may cause a decrement in Ci(Av) which is 
larger than that due to the signal. We use Fi(Av) defined as the 
ratio of the two decrements 



Fi(Av) 



{Cl{0)[l - Kl{Av)]}signal 

{C;(0)[1 - I;(A^)]} Forogrounds 



(34) 



to asses the feasibility of separating the HI signal from the 
foregrounds. This gives an estimate of the accuracy at which the 
Av dependence of the foreground Ci(Av) has to be characterized 

for the signal to be detected. Note we assume that the ^ 
term in eq. (32) can be factored out before considering the decre- 
ment in the foreground. Figure 6 shows the results for the DM 
model and the PR model with R = 16 Mpc. First we note that 
Fi(Av) peaks at the angular scales corresponding to I ~ 10, 000 
(ie. 2 ) and the prospects of separating the signal from the fore- 
grounds are most favorable at these scales. A detection will be pos- 
sible in the range I > 1000, Av < lOKHz and I > 400,A^ < 
lOOKHz for the DM and PR models respectively provided the Av 
dependence of the foregrounds C'i(Av) can be characterized with 
an uncertainty less than order unity. The I and Av range would in- 
crease if the Av dependence of the foreground Ci {Av) were char- 
acterized to 10% accuracy. The largest angular scales ( I < 100 ) 
would require an accuracy better than 1% which would possibly set 
the limit for forthcoming observations. 

The angular modes 1 = 1, 000 and I = 10, 000 correspond to 
baselines with antenna separations of ~ 300 m and ~ 3 km respec- 
tively. This baseline range is quite well covered by the GMRT, and 
also the forthcoming interferometric arrays. This is possibly the op- 
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Figure 5. Angular power spectrum C\ (0) at v = 129 MHz for the two most dominant foreground components, the diffuse galactic synchrotron radiation 
(GS) and the extragalactic point sources (PS) assuming Scut = 0.1 mjy. The expected signal is also shown for the DM model and the PR model with 
R = 16Mpc. 
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Figure 6. This shows Ft (Ai/) (defined in equation 34) for the Aw values shown in the figure. We consider both the DM model(left panel) and the PR 
model(right panel). 



timal range for a detection. A possible detection strategy would be 
to use the behavior of Ci (Av) in the range where Fi ( Au) <g; 1 
to characterize the foreground contribution. This can be extrapo- 
lated to predict the foreground contribution at small Av and any 
excess relative to this prediction can be interpreted as the HI signal. 
A very precise determination of the Av dependence of the fore- 
ground contribution would require a very large Av range in the 
region where Fi(Av) <JC 1, and a bandwidth of ~ 10MHz would 
be appropriate. On the other hand, at I ~ 10, 000 the HI Ci (Av) 
decorrelates within ~ 50KHz [or equivalently Fi(Av) shows a 
considerable drop between 10 KHz and 100 KHz, see Figure 6], 
and it would be desirable to have a frequency resolution better than 
~ lOKHz to optimally differentiate between the signal and the 
foregrounds. A lower resolution of ~ 20 KHz would possibly suf- 
fice at I ~ 1, 000, particularly if the PR model holds. 
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APPENDIX A: CALCULATION OF THE ANGULAR 
POWER SPECTRUM Q (i/i, v 2 ) 

In this section, we present the details of the calculation for the 21 
cm angular power spectrum Ci(vi, v 2 ). The first step would be to 
calculate the spherical harmonic component a ;m of the brightness 
temperature T(v, fi). Using the expression (6) for f/m (k) in equa- 
tion (10), the expression for ai m can be written as 



ai m (v) = T x m J dQYi* m (h) J 



k 



(2tt)3 

A H i(k) + (n-k) 2 A(k)j c - ifcr " (fc - fl) 



(Al) 



which then essentially involves solving angular integrals of the 
forms /dft Y{ m e - ifcr -( fc fl ' and / dfi Y* m (fi ■ k) 2 C - ifcr - (ii fl) 
respectively. Expanding the term e ~ lfc,v ( k ' n ' in terms of spherical 
Bessel functions ji(kr u ), one can show that 



dQ YUn) e-^ *) = 4^H) i j ; (fcr„)F ; ; n (k) (A2) 
Differentiating the above equation with respect to kr v twice 



/ 



dtt (k-n) 2 YT m (n) e" 1 



kr u (k-n) 



= -47r(-i) ; j i "(fcr„)y i * m (k) 
(A3) 

where j" (x) is the second derivative of ji(x) with respect to its 
argument, and can be obtained through the recursion relation 



(2/ + = ^—^- 2 (x) 



+ 



21 - 1 
(/ + !)(/ + 2) 



I 2 



21 - 1 



21 + 3 



21 + 3 

So the final expression of ai m is given by 



31 i x ) 
(A4) 



^Tx m (-i) 1 



d k 



(2tt) 3 " 

x [Am(k)j,(*T^)-A(k)j,"(*T.)] 



(A5) 



The next step is to calculate the the power spectrum 
Ci{v x ,vi) = (ai m (v{) a ; * m (^ 2 )}. The corresponding expression 
is then given by 

Z 2 ) XHl(«l)XHl(«2) 

f d A ki f d A k 2 . 

x ([A H i{zi,\s.i)ji(kir vl ) - A(zi,ki)j,"(fcir l , 1 )] 
x [Am{z2,k 2 )ji(k 2 r V2 ) - A*(z 2 ,k. 2 )j' l '(k 2 r U2 )]) 

(A6) 

where we have put back the redshift-dependence into the expres- 
sions for clarity. Now note that we would mostly be interested 
in cases where ^2 — ZA = Ai/ -C v\. In such cases, one 
can safely assume T(z 2 ) ~ T(zi) and xm(z 2 ) ^ xui(zi). 
Furthermore, the terms involving the ensemble averages of the 
form (AA*) can be approximated as (A(zi, ki)A*(«2, k 2 )) ~ 
(2-7r) 3 (5ij(ki — k 2 )P(2i, fci) and similarly for terms involving 
Am. We can then use the Dirac delta function fo(ki — k 2 ) to 
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compute the k 2 -integral, and thus can write the angular power 
spectrum as 

Ci(Au) = Ci(v, v + Ai/) 

= (4^) 2 f 2 x 2 m J ^y4(ki)y Jm (ki) 

x [jiikr^jiikr^P^ik) 

- {ji(kr„)ji \kr V2 ) +ji{kr„ 2 )j"{kr v )}PA HI {k) 

+ fi'(kr v )j'/(kr„ 2 )P(k)] (A7) 

Using the normalization property of the spherical harmonics 
J dn|F; m (fi)| 2 = 1, one can carry out the angular integrals in the 
above expression, and hence obtain the final result (11) as quoted 
in the main text. 



APPENDIX B: CORRESPONDENCE BETWEEN ALL-SKY 
AND FLAT-SKY POWER SPECTRA 

As discussed in section 2.3, we shall mostly be interested in very 
small angular scales, which corresponds to I 3> 1. For high values 
of I, it is most useful to work in the flat-sky approximation, where a 
small portion of the sky can be approximated by a plane. Then the 
unit vector n towards the direction of observation can be decom- 
posed into fi = m + 6, where m is a vector towards the center of 
the field of view and is a two-dimensional vector in the plane of 
the sky. 

Without loss of generality, let us now consider a small region 
around the pole 9 — * 0. In that case the vector 9 can be treated as 
a Cartesian vector with components {6 cos 4>, 9 sin </>}. This holds 
true for any two-dimensional vector on the sky, in particular U = 
{U cos (f>u, U sin 4>u}- Then the spherical harmonic components 
of T(v, n) [defined in equation (10)] can be written as 

ai m (v) « j a\9Y; m {6,<t>) T(i/,n) (Bl) 

where we have replaced J dtl — * J d9. Now use the expansion 



^(-i) m J m (2^W)e i 



part. This gives a relation between the flat-sky Fourier transform 
T(f, U) and and its the full-sky equivalent ai m {v). 

Using the above relation, we can calculate the power spectrum 



{f(u u U)T> 2 , U')> « \j j^j E (-i) m - m 'e- im *« e im '^' 



X (a2 I r(/,m(vi)02 Ir[ ,/ ira »(f2)) (B6) 

Use the definition (aj m (i>i)aj/ m ; (^2)) = Ci8 H i5 mm > and the 
property 

^e-Mfe-V) =2-k5%\4> u -4> uI ) (B7) 



to obtain 



(f (/a,U)T> 2 ,U')> = 2ttCWK^) &J jf- 5 { V{<t>u-<t>u<) 

(B8) 

The last step involves writing the right hand side of the above equa- 
tion in terms of the two-dimensional Dirac delta function, which 
follows from the expansion 



<^ 2) (U-U') = J 



d0e 



-2iri(U-U')-e 



(B9) 



The exponentials can be written in terms of the spherical harmonics 
using equation (B4): 



4 2 >(u -u')« J do sJj^Y, H)"' 

mm' 

x Y£ Uim (0, tt>)Y^ w , m , (9, 0) e - imfe e im '^' 

(BIO) 

Finally use the orthonormality property of spherical harmonics 
/ d9Y l * m (9, (j))Y Vm , (9, 0) = 5 w 8 mm , and the relation (B7) to ob- 
tain 

<4 2) (U-U') =2n^6 D {4>u-4>u>) (BID 

Putting the above relation into (B8), we obtain equation (20) used 
in the final text. 



im{<t>u—<l>) 



(B2) 



where J m (x) is the ordinary Bessel function. Further, we use the 
approximation for spherical harmonics 



Ylm(e,<t>)^ Jm(W)\U- e™ 
V 27T 



(B3) 



to write 



jj J2(- [ r Y ^u, m (0,^ im * u (B4) 



Then the two-dimensional Fourier transform of the brightness tem- 
perature [defined in equation (16)] will be 

T>,U) = J d6» C - 27riu e T{y, ft) 

* E(- i ) m ^ im ^ / A0 Y ^u, m (0, 4>) T{v, ft) 

= J2(-i) m e~ im * u a2*u, m (v) (B5) 



where we have used the expression (Bl) for ai m in the last 



